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A computational fluid dynamic model was used to simulate the thermal destratification in 
an upright self-pressurized cryostat approximately half-filled with liquid nitrogen and 
subjected to forced sinusoidal lateral shaking. A full three-dimensional computational grid 
was used to model the tank dynamics, fluid flow and thermodynamics using the ANSYS Fluent 
code. A non-inertial grid was used which required the addition of momentum and energy 
source terms to account for the inertial forces, energy transfer and wall reaction forces 
produced by the shaken tank. The kinetics-based Schrage mass transfer model provided the 
interfacial mass transfer due to evaporation and condensation at the sloshing interface. The 
dynamic behavior of the sloshing interface, its amplitude and transition to different wave 
modes, provided insight into the fluid process at the interface. The tank pressure evolution 
and temperature profiles compared relatively well with the shaken cryostat experimental test 
data provided by the Centre National D’ Etudes Spatiales. 


Nomenclature 

a = acceleration 

E = energy 

F = surface force 

k = thermal conductivity 

171 = mass flow rate 

P = pressure 

S = source term 

t = time 

T = Temperature 

u = velocity 

x = lateral displacement 

X = displacement amplitude of oscillation 


Greek 

a = Volume of Fluid fraction 

y = curvature 

p = effective viscosity 

p = density 

<7 = accommodation coefficient 

t = time constant (ramp up of shaking) 

co = shaking frequency 


Subscripts 

eff = effective 

i = interface 

v = saturated vapor 
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I. Introduction 

F uel or propellant sloshing in launch vehicles can have challenging and detrimental effects on the performance of 
the vehicle. Any sudden propulsive maneuvers can instantly induce large motion sloshing behavior. Due to the 
large proportion of propellant mass these vehicles carry, the sloshing of the liquid fuel against the tank walls and 
inertial response of the liquid mass can produce large or critical dynamic loads on the vehicle. In addition, the liquid 
and vapor phase can undergo significant thermodynamic changes that can also produce strong pressure waves and 
variations. Therefore understanding the real-time dynamics and thermodynamics of cryogenic tank sloshing is vital to 
predicting and controlling the tank temperature and pressure response on launch vehicles, or future propellant depots. 

A partially filled tank is most susceptible to sloshing when undergoing lateral vibrational or oscillatory motion. 
The ullage region in cryogenic propellant tanks will typically have significant temperature variation between the tank 
walls and liquid/ullage interface (thermal stratification). Motion of the liquid/ullage interface of sufficient amplitude 
during sloshing events can lead to thermal destratification and a significant decrease in tank pressure. The dynamics 
of the sloshing phenomenon has been studied quite extensively since the 1960’s (see e.g. Refs 1 and 2). Analytical 
studies of sloshing in cylindrical tanks have also been conducted 3 . However the interest in thermal destratification is 
more recent. Thermal de-stratification and pressure reduction from lateral sloshing has been investigated 
experimentally by Lacapere et al. 4 , Arndt et al. 5 and Himeno et al 6 . 

This paper will describe the full three-dimensional CFD simulation and comparison with experimental data of an 
upright 182 mm diameter and 800 mm tall cylindrical tank approximately half filled with liquid nitrogen and subjected 
to forced sinusoidal horizontal shaking. The liquid nitrogen sloshing experiments 4 were performed in 2002 by the 
Advanced Technical Division of Air Liquid as part of the French-German COMPERE (COMPortement des ERgols 
dans les REservoirs) program. Information on the experiment was provided to NASA by Centre National D ’Etudes 
Spatiales (CNES) of France under an Implementing Arrangement on CFD benchmarking for cryogenic propellant 
tank analysis. Improving the capability of computational tools to predict fluid dynamic and thermodynamic behavior 
in cryogenic propellant tanks under settled and unsettled condition is part of the work supported under the previous 
NASA Cryogenic Propellant Storage and Transfer project and under the current NASA Evolvable Cryogenics project. 


II. Computational Model 

The computational (CFD) model solves for the equations of continuity, momentum, and energy. The Volume of 
Fluid (VOF) method is employed to track the liquid and gas phases and their interface in the computational domain. 
The energy and temperature are defined as mass average scalars, e.g.: 
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Here a, p, E are the phase fraction, density, and energy respectively. The phase properties in a two-phase system are 
determined by their volume-fraction averages quantities for density, effective viscosity and thermal conductivity as: 
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The mass conservation of the VOF phase is provided by 
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where S is a source term to account for mass generation and sink sources. 

A Continuum Surface Force (CSF) model proposed by Brackbill et al. 7 was implemented to incorporate the effects 
of surface tension attributed to the VOF phases, and defined by: 
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where y = V • h is the curvature on each of the surfaces (j) on the cell and h is the unit normal vector of the 
surface. 


The ullage was represented by a compressible ideal gas. Simulations were performed treating the flow as laminar or 
turbulent. For turbulent simulations, the k-cd SST turbulence model of Menter et al. 8 with a variable turbulent damping 
value of 10 was used. A laminar simulation was also performed for comparison. The thermodynamic properties of the 
phases was kept constant during the simulation. A second order upwind scheme was used for discretization of the 
turbulence and momentum equations. 


A. Computational Grid 


A full three-dimensional (360 degree) grid was used in order to account for any asymmetric effects, in particular from 
the interaction of the highly unsteady flow in the two phases, swirling surface waves, and the unsteady sloshing 
dynamics. The grid consisted of structured hexagonal cells with some skewness of the cells near the cylindrical wall. 
Three different resolution and grid density configurations were investigated to study their effects on the 
thermodynamic variables and the flow evolution. Figure 1 shows the different computational grids. Here the cell 
dimension refers to the height of the cell which was the controlling dimension specified in the gridding process (the 
arc length of the outer cells on the tank wall were also similarly specified). The other dimensions were automatically 
determined by the gridding software algorithm and were typically smaller than the controlling dimension. The lowest 
resolution grid (left), having 579423 cells, was composed of 3mm cells in the mid-region of the tank to resolve the 
sloshing dynamics and thermal destratification in the interface region. The top and bottom portions of this grid, where 
the flow is expected to be relatively subdued, consisted of 10 mm cells at the very top and bottom rows of cells 
respectively and progressively got smaller as the cells approached the mid-region of finer cells. The medium resolution 
grid (middle), consisting of 741273 cells, was similar to the lowest resolution grid except that the fine 3 mm cells of 
the middle region extended to the top wall, covering the top 3 A of the tank, and the bottom region grid started out with 
20 mm cells on the bottom wall. The highest resolution grid consisting of close to a million cells (984048) had a very 
narrow fine grid region of 1 mm cells in the neighborhood of the free interface and transitions to 3 mm cells on the 
top and bottom regions. 




10 mm cell on 
top with 1 % 
shrinkage on 
consecutive 
cells towards 
mid-section 


I 3 mm cells in 
I mid-section 
r of the 
I cylindrical 
I geometry, 

} Same as on 
top start with 
10 mm cell 
on bottom 



3 mm cells in 
*top 3 /4of the 
cylindrical 
geometry, 


} 20 mm cell on 
bottom with 1 0% 
shrinkage on 
consecutive 
cells 



z 



Refined mid- 
region: 1 mm cell 
1 in center region 
j with 1 .5% cell 
growth to upper 
and lower 
boundaries 

3 mm cells above 
and below the 
refined mid-region 


Figure 1: Three computation grids (coarse to very fine resolution from left to right) 
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B. Non-inertial grid sloshing model 


The tank was modeled in a non-inertial grid and coordinate system. To account for the lateral shaking of the tank, 
source terms were added to the momentum and energy equations that simulated the inertial forces, corresponding 
energy transfer, and reaction forces by the tank. The momentum source term (in the direction of forced lateral shaking) 
is given by*: 


'x-mom = ~ a P ( 5 ) 


The energy source due to forced lateral shaking is given by: 

S = —a pu = S u (6) 

x-energu A' x-mom v ' 


The physical lateral shaking motion of the tank was represented by the product of a sinusoidal and a time- 
dependent function. The lateral shaking ramps up from the initial static position to the maximum amplitude within 
seconds after the shaking is started. Equation 7 gives the displacement in terms of the sinusoidal shaking and the 
exponential ramp up: 


x = X sin (cot) 


-t 

l-e T 


(7) 


X , co, and r are the shaking amplitude, frequency, and exponential time constant respectively. Eq. (7) is incorporated 
into the source terms of Eqs. 5 and 6. The time constant, x, was calculated from the CNES video data of the shaken 
cryostat. 


C. Interfacial mass transfer model 


The kinetics -based Schrage 9 equation is used to model the interfacial mass transfer due to evaporation and 
condensation. The net mass transfer is given as a mass flux vector between phases: 


S„ = m A 
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Here fh is expressed in terms of Schrage’ s relation as: 
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where A- = V<Z , o is the accommodation coefficient, and P and T are the pressure and temperature at the saturated 

vapor state (v) and at the interface (i). This relation incorporates an accommodation coefficient, a, which dictates the 
rate of mass transfer between the phases at conditions near the saturation point. The sign of the accommodation 
coefficient indicates, depending on convention, whether the liquid is evaporating or condensing. 


* Note that source terms are given in terms of the rate of change of the variable of interest. For instance, the source 
term for momentum will have units of mass/length 2 /time 2 . 
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D. Model Setup 

As starting conditions, the tank is initially self-pressurized and thermally stratified by imposing an initial 
temperature field provided by the CNES experimental data. The initial thermal stratification is mostly confined to the 
ullage region of the tank. This initial condition simulates the effects of self-pressurization that a tank may typically 
undergo prior to launch. The computational domain is initially partitioned (patched) so that it is filled with liquid 
nitrogen to a height of 450 mm and gaseous nitrogen above that in the ullage region. The thermodynamic properties 
of both phases are kept constant at saturation values of the initial pressure (2.5 kPa). In addition, several local heat 
loads, shown in Fig. 2, are specified along the cylindrical tank wall that remained constant during the simulation. The 
tank is shaken sinousiodally and laterally at a frequency of 2.1 Hz and a displacement of 3 mm. The 2.1 Hz is close to 
the natural frequency of the partially filled tank system. 4 


adiabatic 



Figure 2: Thermal boundary conditions on the tank. 

III. Results 


A. Pressure evolution 

A set of simulations over the initial shaking period were performed to assess the sensitivity to the effects of the 
physical models and grid resolutions incorporated in the simulations. Figure 3a shows the effects of the interfacial 
mass transfer and inertial shaking. The pressure has been normalized by the initial reference pressure. As seen in 
these plots, the independent contribution from either interfacial mass transfer, i.e. evaporation and condensation, and 
inertial shaking on the pressure was moderate. However, when both sources are combined the effect on the pressure 
is significant. In fact, the decrease in pressure in this case is more than the additive effects of the individual 
contributions, suggesting perhaps a non-linear coupling between these two effects. The effect of grid resolution is 
presented in Fig. 3b which shows a sudden convergence of the higher resolution grids. This seems to indicate that 
the flow and thermal scales were not well resolved when modeled with the coarse grid (which increases to 10 mm 
cells at the top of the tank). 

The effects of the rate of mass transfer at the interface, controlled by the accommodation coefficient in the model, 
was also investigated. The value of the accommodation coefficient produced large effects on the depressurization of 
the tank at low values. Figure 4 shows a large increase in depressurization between the values of 0.0001 and 0.001. 
At values between 0.01 and 0.1, the effects on the tank depressurization were minor. It is clear that the model is 
quite sensitive to small changes in the rate of mass transfer at low ranges of mass transfer rates, and insensitive at 
larger rates. 

The modeling data is compared to the CNES experimental data in Fig. 5. The pressure during the initial phase 
decays more gradually in the experimental case than in the modeling results. Although a time constant (in Eq. 7) 
based on the experimental data was incorporated into the model, the modeling results could not reproduce the more 
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gradual experimental depressurization. Two cases of the modeling data are plotted for comparison. In the case where 
the fluid flow in the liquid and ullage region where modeled as laminar the overall depressurization is slow 
compared to the experimental data. When turbulence is modeled, using a k-cd SST model, the effects on the thermal 
mixing are significant and produces a larger drop in the pressure. In fact the pressure in the turbulent case 
approaches the experimental data at the latter stages, approximately after the 20 second mark. 



(a) 



(b) 


Figure 3: Pressure plots: (a) comparison different models, (b) comparison of grid resolution. 



Figure 4: Pressure time history plots for different values of the accommodation coefficient 
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Figure 5: Pressure plot comparing modeling and experimental data. 

Sloshing dynamics and de-stratification 

One of the first set of observations obtained from the simulation was on the dynamics of the sloshing free surface 
wave which seemed to evolve over different stages. This aspect of the sloshing process has implications on the ullage 
region since the interactions occurring at the interface have a strong bearing on the thermodynamics near the interface 
and throughout the ullage. Three distinct transitional phases in wave structural dynamics are shown in Fig. 6. Here the 
interface is represented by an iso-surface corresponding to a liquid volume of fluid fraction equal to 0.5. The 
computational results incorporates the interfacial mass transfer (with an accommodation coefficient of 0.01) and 
inertial shaking models. The first transition occurred at around 4 seconds after the start of shaking. Secondary wave 
modes were generated along the liquid/ullage interfaces indicated by the folds or creases in the iso-VOF surface. Prior 
to this the interface was smooth with no creases. Slightly later, at 5 seconds, the interface showed signs of wave 
breaking at the periphery where pockets of the liquid nitrogen are ejected from the free surface. Then, at 7.8 seconds, 
a dominant and mainly persistent tangential (swirl) wave mode appears (the grid on the initial free surface is provided 
to show the tilting or precessing of the free surface in the process of swirling). 

1 = 4 sec t = 5 sec t = 7.8 sec 




Secondary Wave Tangential/Circumferential 

wave modes breaking mode 


Figure 6: Sloshing dynamic transitions. 

The temperature contour plots in Fig. 7 show the evolution of the thermal de-stratification. The temperature has 
been normalized by a reference initial liquid nitrogen temperature provided by CNES under the Implementation 
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Agreement. Beginning at the 5 second mark, the initially well thermally stratified ullage starts to show a significant 
convective mixing layer (large scale structures) in the region above the interface. The thermal mixing layer quickly 
grows and propagates upward into the ullage, reaching the top of the tank at around 25 seconds when temperature 
variations there are first seen. At 40 seconds the thermal stratification is significantly altered with maximum 
temperature ratio in the ullage reduced to about 80% of the initial normalized temperature. 

It was also found that the sloshing motion produces different stages of periodicity and intermittency of the temperature 
field. In Fig. 8a the local temperature history at positions within 2 cm of the interface (shown as the 0.4 to 0.47 m 
vertical position values from the tank bottom) show pronounced large amplitude periodic behavior. This periodic 
behavior persists for about 10 seconds slightly after the start of shaking, and then transitions suddenly into an 
equilibrium state characterized by low amplitude random fluctuations. In the central region of the ullage however, at 
a vertical position of 0.6 m, the temperature is characterized by an initial quiescent period of about 6 seconds followed 
by one dominated by a very large amplitude intermittency. After this center-ullage temperature transient, the 
temperature transitions to one exhibiting a lower amplitude fluctuating random behavior. It appears that the thermal 
convective motion evolves on two distinct spatial and time scales. In the mixing layer above the free surface, the 
convective thermal transport evolves more rapidly but is confined to a narrow layer, i.e. smaller temporal and spatial 
scales. By comparison, in the central ullage region the larger spatial convective structures evolve over a corresponding 
longer time scale. 

Figure 8b shows the vertical temperature profiles obtained at different instances of time. It clearly shows that the 
dominant thermal stratification is confined mainly to the ullage region. At the start of the simulation the temperature 
is almost constant in the liquid phase and then becomes slightly stratified a later times. In the ullage, the opposite 
effect, but a larger scale, takes place. It shows more extensive destratification that progresses monotonically with time. 


IV. Conclusions 

The modeling of a partially filled sloshing cryostat tank undergoing thermal destratification was investigated. The 
CFD simulations incorporated models for inertial shaking and mass transfer at the interface. The system pressure 
decayed quickly and at a significant level in both the computational and experimental cases. The turbulent 
computational data with an accommodation coefficient of 0.01 produced the closest fit to the CNES pressure data. 
The sloshing free surface was found to initially oscillate to a single mode and then transitioned into secondary and 
tangential (swirl) wave modes. The temperature field in the ullage was characterized by the propagation of an 
energetic thermal convective layer generated at the sloshing free surface into the upper regions of the ullage. The 
temperature data also showed that sloshing free surface induced two distinct spatial and temporal thermal convective 
scales and periods of intermittency and periodic behavior before reaching an almost equilibrium state. At the end of 
40 seconds the thermal stratification is significantly altered with maximum temperature ratio in the ullage reduced to 
about 80% of the initial normalized temperature. 
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Figure 7: Time evolving temperature contours on the mid plane of the tank. 
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(a) 


(b) 


Figure 8: Temperature plots: (a) time history (b) vertical temperature profiles. 
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